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We review the main steps of the differential equation approach to the analytic evaluation of Feynman graphs, 
showing at the same time its application to the 3-loop sunrise graph in a particular kinematical configuration. 


• The differential equation approach to the 
analytic evaluation of Feynman graph inte¬ 
grals applies to loop integrals defined in the 
by now customary regularization in d con¬ 
tinuous dimensions. We will recall the main 
steps of the approach, describing as an ex¬ 
ample its application to the 3-loop sunrise 
graph in the kinematical configuration of Fig.l. 



Fig. 1. The considered 3-loop sunrise graph. 

• As the very first step, for any given graph one 
has to take all the scalar products which can be 
formed with the external momenta and all the 
loop momenta. In the case of Fig.l there is a 
single external momentum p and three loop mo¬ 
menta ki, i = 1, 2,3, so that, irrespectively of the 
actual structure of the graph, the scalar products 
are the three {p ■ ki), the three kf and the three 
{ki ■ kj) for i 7 ^ j, for a total of nine. One can 
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choose to express some of the scalar products as 
linear combination of the denominators (in gen¬ 
eral, there are several possible choices, equivalent 
for the following); the surviving scalar products 
are said irreducible (perhaps better not trivially 
reducible). One then considers the whole fam¬ 
ily of all the integrals obtained by taking as in¬ 
tegrand all the possible combinations of the P 
denominators of the scalar propagators raised to 
integer powers a,, q = 1,P and of the S irre¬ 
ducible scalar products raised to integer powers 
j3r, r = 1, ...,5”, the continuous dimensional reg¬ 
ularization ensuring the absence of convergence 
problems. The corresponding integrals are func¬ 
tions of the continuous dimension d, of the Man¬ 
delstam variable p^ and of the set of (P -I- S) in¬ 
dices aq,j3r. (For Fig. 1, P = 4 and S = 5). 
For the following discussion, let us further define 
A = J^q(ag - 1), with aq > 1 and B = 
with /3r > 0. 

• As next step, we look for relations between the 
previous integrals. The relations are best pro¬ 
vided by means of the by now famous integration 
by parts identities {ibp-id’s) [1] , obtained by re¬ 
placing in any of the above integrals the integrand 
by all the combinations of the kind {d/dki^^)v^, 
times the same integrand, where i = I, ..,3 and 
Vfj, is any of the external or loop momenta. (For 
the graph of Fig, 1 there are 3 x 4 = 12 different 
combinations {d/dki^^)v^.) 

By carrying out explicitly the derivatives, and af¬ 
ter some trivial algebra, one is left with a com- 
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bination of polynomials (depending in general on 
d, on the Mandelstam variables, and on the inter¬ 
nal masses, always with integer coefficients) times 
integrals of the same family. But the whole orig¬ 
inal integral vanishes, as the integrand is by con¬ 
struction a divergence, so that the combination 
of integrals vanishes as well - an integration by 
parts identity has been obtained. If the original 
integrand belongs to the class {A, B), i.e. if the 
sums of the powers of denominators and numer¬ 
ators are A + P and B respectively, the involved 
integrals belong to the classes {A±1,B ±1). In¬ 
tegrals in which one of the propagator is missing 
{i.e. its exponent is 0) are said to correspond to 
subtopologies, and can be considered as already 
known in a systematic bottom up approach. (It 
is to be observed that other identities of similar 
nature can be obtained by exploiting the Lorentz 
Invariance [2] on the external vectors, but only 
when there are at least three independent exter¬ 
nal vectors.) 

• At this point, it is natural to try to exploit the 
identities for expressing as many as possible of 
the integrals actually needed in any explicit cal¬ 
culation in terms of as few as possible integrals 
of the set - the so-called master integrals. The 
original ideas for building algorithms applying to 
the most general case goes back to Laporta ([2]- 
[5]): the integrals are ordered by giving them a 
weight (that can be done almost at will by giv¬ 
ing a higher weight to the integrals considered 
more complicated; typically, the fully scalar in¬ 
tegral is given the lowest weight), and then a 
“sufficient number” of explicit equations is writ¬ 
ten for them. In the case of Fig. 1, starting 
from the fully scalar integral, corresponding to 
{A = 0, B = 0), one obtains 12 ibp-id’s involv¬ 
ing, besides the original integral, the 4-1-5-1-20 = 
29 integrals with {A = 1, B = 0), (A = 0, B = 1) 
and (A =1,5 = 1) -an apparently runaway sit¬ 
uation, in which the involved integrals outnum¬ 
ber the equations. But that is not the case; an 
elementary combinatorial calculation shows that 
by writing all the ibp-id’s starting from A = 0,1 
and 5 = 0,1 one generates 12 x 30 = 360 equa¬ 
tions, involving integrals with A = 0,1,2 and 
5 = 0,1, 2 at most; as the number of all the inte¬ 
grals with A, 5 = 0,1, 2 is just 336 one has more 


equations than integrals, i.e. an apparently over¬ 
constrained system. That is not the case either 
- simply, the equations are not all linearly inde¬ 
pendent. 

The system can be solved by Gauss substitution, 
eliminating first the integrals of higher weight - 
the method is elementary, even if, given the size 
of the system, it is convenient to use an autom¬ 
atized computer procedure. When that is done, 
it is found that all the integrals appearing when 
dealing with the graph of Fig.l. can be expressed 
in terms of 4 Master Integrals (Mi’s). 

• Let F{d,p^) be any of the Mi’s. Obviously, 

P‘-^e{d.P) = 2p„ 

Substitute in the r.h.s. for F{d,p^) its original 
definition as integral on the loop momenta of a 
suitable integrand. By carrying out explicitly the 
derivative d/dp^ and then contracting with p^^ 
one obtains a combination of integrals associated 
to the considered graphs, which can then be ex¬ 
pressed, thanks to the solution of the ibp-id’s dis¬ 
cussed in the previous paragraph, as combination 
of the Mi’s. From the above equation one has 
therefore expressed the p^-derivative of the par¬ 
ticular MI F{d,p^) in terms of the Mi’s of the 
problem. By repeating the procedure for all the 
Mi’s, one obtains a system of linear first order 
differential equations for the Mi’s. The system is 
in general non homogeneous, as the ibp-id’s can 
introduce integrals where some of the propagators 
are missing. The non-homogenous part involves 
simpler graphs, and can therefore be considered 
as known in a systematic bottom up approach. 

If the Mi’s depend on several independent mo¬ 
menta and correspondingly by several Mandel¬ 
stam variables, by acting on the Mi’s with all the 
combinations and minor algebraic rear¬ 

rangements one can obtain a system of differential 
equations in any of the Mandelstam variables [6]. 
It is to be recalled here that the solutions of the 
ibp-id’s can be almost immediately rewritten as 
differential equations in the internal masses [7]. 
Finally, a linear system of differential equations 
can always be rewritten as a single differential 
equation of suitable higher order for any of its 
unknow functions. The coefficients appearing in 
the equations are in any case polynomials (or ra¬ 
tional factors, when dividing the equations by the 
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polynomial multiplying the highest derivative). 

• In the case of the special kinematical situation 
depicted in Fig.l, one of the internal masses van¬ 
ishes, two internal masses take the same value 
M, the fourth the value m and the Mandelstam 
variable is timelike and takes the value i.e. 

= — rn^. We rescale all the momenta by M, 
further put m = Mx and define 

= (Z)3(S2) j ^‘'^3 X 

^l(^2 + l)(fc3 + 1)[(P — fcl — ^2 - + X^] ’ 

where C{d) = r(3 — d/2) is a normal¬ 

ization factor, introduced for convenience, with 
the limiting value C(4) = 1 at d = 4. The above 
discussed equations in p'^ and the masses give for 
^{d,x) the following 3rd O.D.E., exact in d: 

3 

+Void,x) $(d,x) = 

N{d,x), (2) 

where 


the first equation involves only 
and is to be solved for independently of 

the other coefficients; the second equation is an 
equation for involving in the 

non-homogeneous term, and so on. The general 
equation reads 

^ Po{x) 

2 ,i 

= + Qoix) 

1^1 

+ \Riix)^ + Roix)]^^^-^\x) 

L ax J 

+So{x)jj^\n(>^+^\x) (5) 

where = 0 if /c < —3, and Pi{x), Qi{x), 

Ri{x), S'o(a;) are simple polynomials in x whose 
explicit expressions can be easily derived from 
Eq.(2). 

All the above equations have the same form 
V{x)^^'^\x) = [x) , ( 6 ) 


'P 3 {d,x) = x^{l - x){l + x) ] 

V 2 {d,x) = [2-(d-4)]a;-b [2-k5(d-4)]a;3 ; 

Pi{d,x) = -[6-b8(d-4)-k2(d-4)2] 

-b[2-4(d-4) -6(d-4)2]x2 ; 
Po{d,x) = - [8-I-14(d — 4)-I-6(d — 4)^]a; ; 

1 

W-) = (jvyi—^ (=>» 

• The direct solution of the differential equation, 
for arbirary d, seems out of reach. But as, in 
any case, one is interested in the d ^ 4 limit of 
the solution, one can expand $(d, x) in Laurent 
series in (d — 4) and rewrite the equation in terms 
of the coefficients of the expansion. Let us take as 
known, for simplicity, that in our case the leading 
(most singular) term is of order l/(d—4)^, so that 
the expansion of 4>(d, x) reads 

$(d,x) = (d-4)'=$('=)(x) . (4) 

fc>-3 

By correspondingly expanding Eq.(2) in (d — 4), 
one obtains a system of chained equations for the 


where X>(x) is the differential operator applied to 
$(*)(x) in the l.h.s. of Eq.(5) and N^'^^{x) stands 
for the non-homogeneous part, to be considered 
as known when dealing with the order k in the 
(d — 4) expansion. The homogeneous part, which 
is the same for any k, can be written as 



• We will try to solve the equations for the coef¬ 
ficients of the (d — 4) expansion by means of the 
method of the variation of constants by Euler. 
To that aim, we need the solutions of the homo¬ 
geneous equation, Eq.(7). 

Up to here, our approach applies to virtually any 
Feynman graph integral; for continuing one has 
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to solve the homogeneous differential equation in 
d = 4 dimensions. Algorithms for the solution 
of differential equations in the general case are 
unfortunately not available. But as a more opti¬ 
mistic remark, let us point out that in almost all 
the cases considered sofar (see for instance [8], 
[9], [10]) the homogeneous equations came out 
to be rather simple, with solutions which could 
be expressed in terms of almost elementary func¬ 
tions. The same happens also in the present case. 
Indeed, it is immediately seen that a first solution 
of Eq.(7) is 

(f)i{x) = (1 ; 


by putting 4>(x) = (j>i{x)'ip{x) one finds a second 
order equation for ip'ix) having the simple solu¬ 
tion 

= {l-x^f/x^ , 

and finally, putting ^ one finds 

at once 

^'{x) = -xV{l-x^f . 

The corresponding solutions are 


02 (x) = 

1 (1 — x^)(l — x'^) 

~ 2 ^2 


2{l - x^)H{0-x) , 


03 (x) = 

3 (1 — x^)(l — x^) 


[id(-l;x) + id(l;x)] 


512 a;2 

+ [id(0,-l;x) + id(0,l;x)] 

1 (x^ -I- 3)(3x^ -I-1) 

~ ^ X ’ 

where the functions H are Harmonic Polyloga¬ 
rithms or HPLs [11]. Referring to [11] for more 
details, let us recall that the HPLs depend on 
an argument, say x, and on a vector of indices, 
say b, whose components take any of the values 
(1, 0, —1), and whose number is called the weight 
(w). At re = 1, they are defined as 


H{0;x) = ln(x) , 

r dx' 

H{l;x) = -\n{l-x)= ; 

Jo ^ ^ 

r dx' 

H{-l-x) = Hl+x)= 

Jo -L “r 3; 


at higher weight, ru > 1, if the vector of the in¬ 
dices is written as (a, 5), where 5 is a vector of 
w — 1 components, they fulfill 

■^H{a,b;x) = f{a,x)H{b]x) , (8) 

with 

/(l,x) = 1/(1-x), 

/(0,x) = 1/x , 

/(-l,x) = 1/(1-fx). (9) 

The Wronskian of the 0i (x), as expected by direct 
inspection of Eq.(7), is 

H^(x) = —(1 — x)^(l -I- x)^/x^ . 


In terms of the solutions of the homogeneous 
equation, the general solution of the inhomoge¬ 
neous equation is given by 



dy 


Cafcc ^bc{y') 

W{y) 


N^'^\y) 


Mbc{y) = <pb{y)Jciy) - Jbiy)My), 


( 10 ) 


where {a,5, c} is a permutation of {1,2,3} and 
= 1,2,3) are three integration constants 
to be fixed, order by order in k, to match the 
boundary conditions identifying <I)((i, cc) Eq.( 1). 
• Let us look at Eq.(lO) for the first value of k, 
which is fc = —3. From Eqs. (5) one gets 




1 (2 _ 1_\ 

2\a; 1 — X l-|-xy 


when it is inserted in Eq.(lO), after partial frac¬ 
tioning and some trivial integration by parts, the 
remaining non trivial integrals involve only one 
of the three factors 1/x, 1/(1 — x), 1/(1 -I- x), the 
three factors /(a, x) of Eq.(9) and polylogarithms 
of various weight: but those (indefinite) integrals 
can be carried out at once by using Eq.(8). It 
follows that N^~‘^'>{x), which contains 4>*^“^)(x), 
is also a combination of powers of the monomials 
X, (1 — x), (I -I- x) and polylogarithms, so that the 
procedure can be iterated at will up to any de¬ 
sired value of k, the result being, at all orders, a 
combination of harmonic poly logarithms. 
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• To complete the job, we must fix the sofar arbi- 

(k) 

trary constants <i>a , a = 1,2,3. To that aim, let 
us observe that the behaviour of the most general 
solution of Eq.(2) is of the kind 

4 

x) = ^ (l + + ....) , 

where the four possible values of the exponents 
are ai = 0, a 2 = —{d — 2), = {2d — 5), and 

a 4 = {d — 2), as can be easily verified by sub¬ 
stituting the above expansion in Eq.(2). On the 
other hand, by direct inspection of x) Eq.(l), 
one finds that the behaviours a 2 and are ruled 
out, as the integral does not diverge in x in the 
a; ^ 0 limit when d > 2; this implies that terms 
like 1/x^ or 1/x (or more in general terms with 
odd powers of x) cannot be present in the expan¬ 
sion around x = 0 in the d —> 4 limit. Similarly, 
one finds by inspection that 4>(d, x) is analytic in 
X for X —> 1; this implies that terms in ln(l — x) 
cannot be present in the expansion around x = 1 
in the d —> 4 limit. Those conditions are sufficient 
to fully determine the integration constants. In¬ 
deed, we evaluated explicitly Eq.(lO) up to fc = 5, 
and by imposing the proper behaviours for x —> 0, 
X ^ 1 we fixed the constants $ 3 *^ for 

—3 < fc < 3, (d) 3 ^^ being 0 for any k) and $ 2 ^^ 
The constants $ 3 ^^ and $ 2 ^^ remained still 
undetermined - to fix them, one has to impose 
the boundary conditions to further terms (x) 
with k = 6,7. 

The explicit analytic solution up to fc = 3 in¬ 
cluded, which we obtained in that way, involves 
HPL’s of argument x and up to weight w = 6. 
The resulting expression is unfortunately too long 
to be listed here. More details can be found 
in [10]. 

• Having the full dependence on x in closed ana¬ 
lytic form it is immediate to obtain the values at 
X = 0 and x = 1. The x = 0 values can also be 
obtained by direct integration of Eq.(l); they can 
be used as a check of the calculation. 

The values at x = 1 are more interesting. They 
correspond to the quantity Jn of [3], where they 
were given up to order (d —4)^; from our solution 
(and from the table [12] of the values of the HPL’s 
at X = 1), we could evaluate them up to order 


(d — 4)® included; in fact, given the structure of 
the solutions of the homogeneous equation, the 
still undetermined constants are not needed for 
the X = 1 value. The results up to (d — 4)^ are in 
full agreement with [3]; the new terms, namely 
the coefficients of order (d—4)"* and (d—4)® match 
the numerical values of the quantity I 12 of [13]. 
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